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We have adapted classical molecular dynamics to study 
the structural and dynamical properties of amorphous silica 
surfaces. Concerning the structure, the density profile ex- 
hibits oscillations perpendicularly to the surface as observed 
in liquid metal surfaces and the pair correlation functions as 
well as the angle distributions show features (absent in the 
interior of the films) that can be attributed to the presence of 
2-fold rings which are perpendicular to the surface. From the 
mean-squared displacement of the non bridging oxygen atoms 
we find that in the interior region they move perpendicular 
to the surface while they move parallel to it in the surface 
region. 

PACS numbers: 61.43.Fs, 61.20.Ja, 68.35.Bs 



I. INTRODUCTION 

Most studies on surfaces have been focused on their 
electronic aspect. The goal was usually to analyze mate- 
rials used in electronic devices which characteristics de- 
pend on the surface properties. As a consequence, crys- 
talline surfaces are currently well known and a lot of at- 
tention has been devoted to the doping and transport 
properties of these surfaces as well as to the surface recon- 
struction. In the last years, many studies, theoretical 0] 
or experimental Q have also been devoted to the analysis 
of liquid metal surfaces. In particular, the existence of an 
atomic layering in the density profile has been found. A 
great interest in understanding amorphous surfaces exists 
1^, because of their participation in the glass formation, 
and therefore several numerical studies are now available 
on this subject [^^. Pioneering work in using classical 
potentials to describe amorphous silica surfaces is due to 
Garofalini et al. who studied mainly the presence of de- 
fects (non bridging oxygens (NBOs), small-sized rings..) 
at the surface Q . These defect sites are the most chemi- 
cally reactive ones but experimental studies have, so far, 
not clearly shown their existence [§,0. An alternative 
way to learn more about the surface is to study the chem- 
ical reactions taking place on it ^ and in particular the 
interaction between water and amorphous Si02 surfaces 
has been the topic of both experimental P,[To| and nu- 
merical studies ||ll|,|l2| . In most of the numerical studies 
available in the literature, the total number of atoms did 
not exceed 800 atoms and only a fraction of these atoms 
were considered as "mobile" (i.e. participating in the 
formation of the surface). (It should be noted that very 



recently (in fact during the process of writing the present 
article) a molecular-dynamics study of silica surfaces has 
been proposed involving several thousands of atoms [^). 
Moreover the long-range Coulomb forces were either cut- 
off 1^ or treated in a pseudo 2D geometry |13 , not taking 
into account exactly the reduced dimensionality of free 
silica surfaces. A recent study of the influence of long- 
range forces on the structure and dynamics of a model 
silica glass has shown that over-damped Coulomb inter- 
actions can alter the quality of the results . 

In the present work the molecular dynamics force cal- 
culation scheme has been adapted to the specific 2D 
geometry of samples with free surfaces. Our aim was 
mainly to study in detail the structural properties of 
amorphous silica surfaces with respect to what is already 
known for the corresponding bulk samples and to test 
the quality of the widely used "BKS" potential proposed 
by van Beest et al. This classical pair potential de- 
scribes quite well the structural and vibrational 
JlSf , as well as relaxational [|l9| and thermal |2^] proper- 
ties of both, supercooled viscous liquid and glassy bulk 
silica samples: it is therefore justified to address its ef- 
fectiveness in the case of "bi-dimensional" silica systems 
even more so since in a recent study based on ab initio 
simulations, Ceresoli et al. have put into question the 
structural description of the surfaces obtained with the 
BKS potential ||2l| . In addition to the work already pro- 
posed by Roder et al. Q who used the same potential 
but aimed their study mainly on silica clusters, we show 
that the effect of the surface is firstly refiected in the 
density profile which exhibits stratification, explained as 
a tendency to have more short range order at the sur- 
face. We then analyze in detail the ring size distribution 
as well as the orientation of the rings with respect to 
the surface, and show that the number of two-fold rings 
(whose signature is visible in the radial pair distribution 
functions and the angle distributions) is coherent with 
the experimental estimates and that they are positioned 
perpendicularly to the surface. By analyzing the influ- 
ence of the annealing time on this effect, we show that it 
is a true surface reconstruction and not an artifact due 
to out-of equilibrium phenomena. Finally, in order to 
elucidate the diffusion process at the surface, we present 
the mean square displacement of the bridging and non- 
bridging oxygen atoms. We find no signature of a liquid 
like behavior at the surface but the dynamics of the NBOs 
seems to be anisotropic. 
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II. SIMULATIONS 



with 



We have performed classical molecular dynamics sim- 
ulations starting with cubic samples of edge length 
L=35.80 A, containing 3000 particles. These samples 
were used previously to study the properties of bulk silica 
therefore their mass density is p ~ 2.18 g/cm'^, which is 
very close to the experimental value of 2.20 g/cm^ [ p2[ , 
and the periodic boundary conditions were assumed in 
the three space directions. The potential used to de- 
scribe the interactions between the atoms is the two-body 
"BKS" potential. After the use of this potential to study 
densified silica samples |2^, we apply it here to study 
silica surfaces assuming this potential is still valid in the 
case of Si02 free surfaces. Of course this description can 
never be as accurate as an ab initio study, nevertheless 
we can simulate much larger samples and typical struc- 
tural features should not depend on the details of the po- 
tential. The classical equations of motion are resolved by 
using the Velocity- Verlet algorithm with a timestep At = 
0.7 fs. The initial configuration is given by a bulk sam- 
ple of amorphous silica at T = K. This configuration 
has been obtained by quenching a well equilibrated liquid 
sample around 7000 K with a quench rate of 2.3 x lO^'' 
K/s. The free surfaces are created by breaking the peri- 
odic boundary conditions along the z-direction, normal 
to the surface, thus creating two free surfaces located 
at L/2 and —L/2. Doing so, the Ewald summation for 
the calculation of the long range forces has to be mod- 
ified to take into account the loss of periodicity in the 
z-dircction. There exists no standard method since sev- 
eral possibilities based on either a modified Ewald sum- 
mation (l^j2^ or multipole expansions are proposed 
in the literature. A purely two-dimensional method ex- 
ists which nevertheless necessitates the adjunction 
of charged plates. Here we have chosen a strictly two- 
dimensional approximate technique which has the advan- 
tage to be very simple In the three dimensional 
Ewald summation technique, one considers all the point 
charges surrounded by a spherical distribution of charges 
that has same magnitude but opposite sign Then, 
the potential is separated into two terms. The first term 
Vc corresponds to the potential due to the point charge 
distribution minus that of the screening spherical charge 
distribution and is quite rapidly convergent in real space. 
The second one Vy is the potential created by the spher- 
ical distribution and is rapidly convergent in the recip- 
rocal space. Since the real space term is a short range 
potential with a spherical symmetry, there is no need to 
change it for the two-dimensional summation. But, of 
course the reciprocal space potential has to be modified. 
While for the 3d Ewald summation the expression for the 
reciprocal space term of the potential is: 
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In our approximate two-dimensional method it becomes: 
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Here A: is a two-dimensional vector and the summation 
runs over the two integers Ux and Uy. Although formula 
(3) looks like a straightforward two-dimensional tran- 
scription of formula (1) it is only approximate. The ex- 
pression would have been exact if, instead of a Coulomb 
potential with a spherical symmetry one would have con- 
sidered a cylindrical potential. In particular, one prac- 
tical consequence is that the k-contribution to the z- 
component of the force is zero, since the scalar product 
k.fi is independent of the third coordinate Zi of fl. We 
have used the characteristic constant k = 5.0/L, where 
L is the cubic box size, and considered 49 k-vectors in 
reciprocal space to insure a relative error smaller than 
2 X 10~^ for the potential energy. The removal of the 
periodic boundary conditions in the z-direction breaks 
the bonds between the atoms at the "bottom" and those 
at the "top" of the film. This will lead to bond rear- 
rangements near each free surface and of course increase 
the temperature of the system. It is also worth notic- 
ing that we let the position of the atoms adjust freely in 
the z-direction and therefore we are not stricto sensu in 
a microcanonical ensemble. To justify a posteriori the 
approximations and check the results we compared this 
method with the one in which the box length in the z- 
direction, , is artificially increased in order to simulate 
a pseudo 2D system. We find that the two techniques 
give similar results when is five times larger than the 
box length in the other 2 directions which is in contrast 
with the study of Bakaev et al. [p^ in which was 
taken only twice as large as in the x and y directions. In 
any case the method used here is computationally much 
more efficient since the number of k-points included in 
the Ewald summation has not to be increased. It is also 
faster than the more frequently used technique proposed 
in Ref. jlj] which is not a strictly 2D technique and which 
is computationally very costly. 

The heating of the slab mentioned earlier will tend to 
evaporate some oxygen atoms located at the surface, and 
to avoid this vaporization, it is necessary to control the 
temperature of the slab. This has been done at lOOOK 
by rescaling the velocities during 30000 time steps. This 
temperature is below the simulated glass transition tem- 
perature as will be shown in a forthcoming study. Then 
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the system was allowed to relax 30000 supplemental steps 
and finally we performed simulations during 60000 sup- 
plemental steps to collect the results, the system tem- 
perature remaining almost constant around 1000 K in 
this time range. Such a simulation represents 45 days of 
computer time on the latest IBM-SP2 processor and in 
order to improve the statistics of the data we averaged 
our results over ten independent liquid samples. 



III. RESULTS 

In figure 1 we have reported the averaged density pro- 
file of our samples in the z-direction. Each point (filled 
circle) of the figure corresponds to the density p{z) cal- 
culated within a slab parallel to the surface of area Lx L 
and of width Az = 0.4475 A and reported as a function 
of the mean slab position z. We have averaged the results 
over the two slabs located at z and —z. For comparison, 
we have reported the profile in the x-direction calculated 
within a slab of same width Aa; = 0.4475 A but of area 
Lx Lz with Lz — 26.85 A, taken sufficiently smaller than 
L to avoid edge effects due to the presence of the free sur- 
face in the z— direction. In both cases, we get an average 
bulk density of about 2.3 g/cm^, significantly larger than 
the density p = 2.18 g/cm'^ of the original cubic box with 
periodic boundary conditions. This increase of the den- 
sity is due to the relaxation of the system (the system 
contracts along the z direction) when creating the free 
surfaces since it has been shown that the equilibrium den- 
sity of the BKS potential is larger than the experimental 
one |2^ (this means that the initial sample with PBC had 
a negative pressure). This shows that creating free sur- 
faces is a physical way to obtain the equilibrium density 
and thus permits also to validate a posteriori our modus 
operandi. For the profile in the z-direction, one can use 
a hyperbolic tangent fit derived from the Van der Waals 
theory of surface tension |30|. Although this formula is 
generally used to study the thermodynamic properties of 
simple liquids, some characteristics of the density profile 
of our amorphous samples can be determined. This fit 
uses three parameters, po the mass density in the under- 
lying bulk, zq the position of the surface and d the surface 
thickness: 
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In fact, this formula fits very well the vapor- liquid inter- 
faces which exhibit symmetric density profiles [Q. But 
here it does not take into account the asymmetric den- 
sity maximum near the free surfaces (see figure 1) and 
therefore the density profile is not very well fitted near 
this maximum. As a consequence, if the parameter po 
(po = 2.36g/cm^) gives a good quantitative value for 
the bulk density, the parameter d (d = 1.18 A ) gives 
only a rough (underestimated) value of the actual sur- 
face thickness. As observed for liquid surfaces [pTI], we 



have an asymmetric peak with a broad tail in the den- 
sity profile at the surface and a layering behavior. Our 
goal is to identify the participation of the stratification in 
the oscillations observed in the density profile along the 
z-dircction. However when comparing the density pro- 
files in the z and x directions in figure 1, one can hardly 
see a difference in the oscillations around the plateau. 
One way to address this point is to compute the stan- 
dard deviation of these data over different samples. If 
the oscillations are only due to statistical fluctuations, 
the standard deviation a of the density should behave 

— 1/2 

like Ng ' , where Ns is the number of statistical inde- 
pendent samples over which the density profiles are aver- 
aged. Along the x and y axis, we find that the standard 



deviations behave like N] 
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while for the z-direction 



they are quite independent on Ns- Therefore the oscil- 
lations along the z-direction are most likely a signature 
of a layer formation whereas the oscillations observed in 
the x-direction are only due to statistical density fiuc- 
tuations. This stratification is likely to be enhanced at 
the free surfaces and to extend in the underlying bulk. 
Computing the radial pair distribution for atoms within 
a sphere of radius R — L/A, whose center is the center 
of the the initial cubic box, we found that it is similar to 
that observed for bulk amorphous silica samples. This 
is an evidence that the layering behavior is not due to 
crystallization. The inset of Fig. 1, exhibits the contri- 
bution of each type of atom to the mass density. Since 
the stoichiometry is respected, the quantity An — 2 [Si] - 
[O] is oscillating around a zero mean within the underly- 
ing bulk. Near the surface, the apparition of two peaks 
refiects some stoichiometry breaking. The silicon atom 
excess under the surface and the oxygen atom excess at 
the surface are respectively represented by the positive 
and the negative peak. This surface segregation appears 
not to be model dependent since it has been found also 
by Roder et al. as well as by Litton et al. with 
another potential. 

In figure 2 we present the radial pair distribution func- 
tions (RPDFs) for the different pairs of atoms in order 
to investigate the structure of our samples. The space 
has been divided into two regions, the "surface" layer 
which extends 6 A below the z position of the outermost 
atom, and the "interior" which represents the rest of the 
sample. The width of the surface layer has been chosen 
to encompass the broad peak in the density profile visi- 
ble in figure 1 and is coherent with previous studies 
Also shown in figure 2 are the RPDFs obtained in a bulk 
silica sample. As a general trend, the RPDFs in the in- 
terior are very similar to the ones obtained in the bulk. 
Nevertheless differences can be seen especially in the O- 
O RPDF, where a slight shift towards smaller distances 
can be seen. This is a consequence of the relaxation to- 
wards a higher equilibrium density. The major difference 
between the surface and the interior can be seen in the 
Si-Si RPDF, since a small shoulder appears around 2.5 A 
. This shoulder which can be attributed to 2-membered 
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rings as we will see below, is absent in the interior or 
bulk curve since this structural unit is not present in the 
bulk |ljjr^,^ . In figure 3 we show the distribution func- 
tions of the angle between an oxygen (silicon) atom and 
its 2 silicon (oxygen ) neighbors (two atoms arc neigh- 
bors if their separation is smaller than the location of 
the first minimum in the corresponding RPDF) where 
we distinguish again the surface, the interior and the re- 
sults obtained for bulk silica (this is a further test of the 
validity of the BKS potential which does not contain 3- 
body terms). Concerning the 0-Si-O angle distribution, 
there is almost no difference between the interior and the 
bulk. The distribution at the surface is slightly broader 
and shows a second peak around 80° which is absent in 
the bulk. Concerning the Si-O-Si angle distribution we 
firstly note a slight shift of the whole distribution towards 
smaller angles compared to the bulk distribution which 
is again a consequence of the density relaxation. Next 
we note a strong shift of the surface distribution towards 
smaller angles by more than 10° which is coherent with 
the peak in the density profile observed in figure 1 and an 
important second peak around 100° which is absent in the 
bulk and interior curves. The existence of the additional 
peaks in the surface angle distributions is in agreement 
with previous studies [Q,!). The location of these peaks 
is different from the angles 0-Si-O and Si-O-Si expected 
in 2 mcmbcrcd rings |pj|: 94 and 86° respectively. Here 
we find 80 and 100° which indicates a weakness of the 
pair potential on this particular point since in a recent 
ab initio study on silica clusters Lopez et al. found 91 
and 89° respectively Q . 

So far we only have the signature of the presence of 2- 
membered rings at the surface therefore it is justified to 
address the influence of the free surfaces on the variation 
of the ring size distribution. Such quantity is of interest 
since it can be extracted from infrared and Raman spec- 
troscopy. In particular the highly strained 2-membered 
rings result in infrared-active Si-0 stretching modes at 
888 and 908 cm^^ 36|. The ring size distribution has 
been determined using the algorithm described in [ p^ , 
where a ring is defined as the smallest loop starting from 
one oxygen atom Oi nearest neighbor of a given Si atom 
and ending on O2 another of its nearest neighbor. Then 
the size of a ring is given by the number of its consti- 
tutive Si-0 segments. To get the distribution of each 
ring size along the z axis in the 2 regions (interior and 
surface) we compute the probability P„ for a given Si 
atom, whose coordinate along the normal direction to 
the surface is in a given region, to be a member of a 
n-fold ring. These results are reported in the top part 
of figure 4 and compared to the bulk distribution. We 
calculated the ring size distribution in the surface layer 
on the initial 2D samples (not relaxed) and on the 2D 
samples annealed for 42ps (relaxed) in order to see the 
influence of the annealing time on the results. Firstly it 
should be noted that the distributions in the interior re- 
gion (relaxed or not relaxed) are very similar to the bulk 
distribution (which is a broad Gaussian centered around 



six-fold rings [Q) and are therefore not reported in fig- 
ure 4a. The situation is different at the surfaces where 
the global trend after the annealing period is a decrease 
of the probability for a Si atom to be a member of a 
large ring (as expected, since a lot of bonds are broken) 
together with an increase of the probability to be a mem- 
ber of a small ring, more precisely, two-fold, three-fold 
and four-fold rings. A particular case of this behavior is 
the inversion of the proportion of five-fold and six-fold 
rings at the annealed surface. Energetically, there is no 
significant difference between these two ring sizes, there- 
fore if six-fold rings are broken at the surface (because 
they are relatively large) a fraction of these rings will 
become five-fold and five-fold rings will become predom- 
inant: in a sense this corresponds to the reconstruction 
of an amorphous surface. An other significant feature 
is the relatively large proportion of two-fold rings at the 
annealed surface whereas it stays very weak in the bulk 
and at the freshly created surface. More precisely we 
find an average density of 0.25 rings per nm^ which is in 
the range of the experimental values which vary between 
0.1 and 0.4 rings per nm^ pC| , |3^ , |35[ | . This proportion is 
smaller than the one found in refW (0.6 rings per nm^), 
where the statistics were done at 3000K and for silica 
clusters at equilibrium (which is certainly not the case 
here). 

Together with the ring size distribution we have also in- 
vestigated the ring orientations by computing < cos^ 9 > , 
for a given ring size n, within a given region of the sample, 
where 9 is the angle between the normal to the surface 
and to the ring. The results for the surface region are 
reported in the bottom of figure 4 as a function of the 
annealing time, and compared with those obtained in a 
3D sample (which are similar to those in the interior re- 
gion). In the bulk, all the n-fold rings are oriented in an 
isotropic way since the value of < cos^ 6* > is close to i . 
In the freshly created surface region, the 2, 3 and 4-fold 
have still an isotropic orientation, while the larger rings 
are oriented rather parallel to the surface since < cos^ 9 > 
is slightly larger than i . In the annealed surface layer we 
note that the rings larger than the three-fold ones have a 
strong tendency to orient parallel to the surface. On the 
contrary the two-fold rings are now perpendicular to the 
surface since < cos^ 9 > is much smaller than i. This is 
a result of the annealing and can therefore be extrapo- 
lated to a film in equilibrium in which two-fold rings are 
certainly oriented perpendicular to the surface. This is 
coherent with the ab initio study of Ceresoli et al. l2l| ] 
even if they had only one two-fold ring at the surface. 
After the description of the structure let us now turn in 
the final part of this work to the diffusion mechanism in- 
side our thin silica films. The motivation of this study is 
connected to the idea that similarly to cn^stalline faces 
which develop a microscopic liquid film the surface 
layers of an amorphous sample show an enhanced diffu- 
sivity compared to the bulk [ ^2|j3S| ]. The simplest way 
to address this question is to compute the mean-squared 
displacement (MSD) {r^{t)) = {\r,{t)-r,{0)\'^) for agiven 
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particle type i. This is what is represented in figure 5 for 
the oxygen atoms (the behavior of the sihcon atoms is 
quahtatively similar as shown in Ref. [^). The interior 
and surface regions have been determined with respect to 
the initial position of the atoms. In figure 5 we have also 
distinguished the behavior of the non-bridging oxygens 
(NBOs) from that of the bridging oxygens (BOs). At 
the surface we find approximately 15 % of NBOs while 
in the interior only 0.5 % of the oxygen atoms are single- 
coordinated, which is in good agreement with a recent 
MD study of nanoporous silica |40j, especially for the 
number of NBOs in the surface region. We have also cal- 
culated the components of the MSD that are parallel ( 
1.5 X {x^ + y^)) and perpendicular (3 x (z^)) to the sur- 
face, and these components are also represented in figure 
5. In order to study the infiuence of the annealing time, 
the MSD has been calculated for 42ps at the beginning of 
the relaxation process (thin lines) and after an annealing 
time of 105ps (bold lines). For the BO there is no dif- 
ference in the MSD between the interior and the surface 
region. After the annealing period, the displacement of 
the atoms is smaller and more isotropic. Hence for a fully 
equilibrated film the dynamics of the particles inside the 
sample is expected to be almost isotropic. 
For the NBOs the situation is different. Firstly as ex- 
pected, the NBOs diffuse more than the BOs. Next, after 
annealing, the dynamics appear to be more anisotropic 
in both regions (interior and surface). In the interior re- 
gion the perpendicular MSD is more important than the 
parallel one as if the NBOs were attracted towards the 
surface. On the contrary in the surface layer the NBOs 
move parallel to the surface rather than perpendicular to 
it. This can be understood since most of the Si-0 bonds 
point outward of the surface and since it is easier for 
an NBO to have a bond "bending" motion rather than 
a bond "stretching" motion. Therefore concerning the 
NBOs, we find that their dynamics is anisotropic both in 
the interior and the surface region. Since this anisotropy 
is enhanced after annealing one can expect that this be- 
havior is also true in an equilibrated sample. 



IV. CONCLUSION 

We have used classical molecular dynamics to study 
the structural properties and the mean squared displace- 
ment of the particles inside an amorphous silica film. The 
first step in this work was to adapt our computational 
method to samples without periodic conditions in one 
spatial direction. This has been done using an original 
method which is fast and relatively accurate. With this 
technique we have calculated the density profile along the 
normal direction z to the surface. We have found that 
the system is contracting since the equilibrium density 
of the potential is greater than the mass density of the 
initial amorphous silica sample. This relaxation is rather 
fast and creating a surface seems therefore to be a physi- 



cal and efficient way to obtain the equilibrium density of 
a given system. This shows also that the BKS potential 
describes properly the physics inherent to the creation 
of free surfaces. Like in liquid metals, the density profile 
exhibits a layering behavior that seems to be enhanced at 
the surface. This feature was interpreted as a sign of an 
increase of the short range order near the surface. The 
asymmetric peak near the surface is a sign of a stoichiom- 
etry breaking: while the silicon atoms are in excess in the 
inner part of the surface the oxygen atoms are in excess 
in its outer part as already mentioned in previous studies 
. We have calculated also the radial pair distribution 
functions and the angle distributions and we found the 
signature of the presence of 2-fold rings at the surface 
similarly to previous studies Therefore we have also 
focused our attention on the ring size distribution and 
ring orientation along the z-direction. We see that at the 
free surfaces, the most energetically unfavorable configu- 
rations (small rings) are enhanced and correlatively, the 
large rings show a tendency to disappear. A particular 
case of this feature is the inversion of the proportion of 
five-fold and six-fold rings. Also the orientation of the 
small-sized rings with respect to the surface (2-fold rings 
are perpendicular to the surface) agrees well with a re- 
cent ab initio study |^^. Concerning the MSD of the 
oxygen atoms, we find that the dynamics of the BOs is 
unaffected by the presence of the surface. On the con- 
trary the dynamics of the NBOs is clearly anisotropic: 
in the interior (far from the surface) they move rather 
perpendicular to the surface while in the surface region 
they move parallel to it. 
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FIG. 1. Density profile p{z) in the direction z, perpendic- 
ular to the surface, from 2 = (median cut of the simulation 
box) up to the outer part of the surface. The filled circles 
correspond to the density calculated within slabs of width 
Az — 0.4475 A for a sample size of L = 35.8 A. The densities 
for slabs of opposite coordinates z have been averaged. The 
curve represents the fit of p{z) using eq. (5). In the inset 
the quantity An = 2 [Si] - [O] is plotted as a function of z, 
where [Si] and [O] are the numbers of Si and O atoms per 
unit volume calculated within the slabs. 
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FIG. 2. Radial pair distribution functions in a bulk silica 
sample (thin solid), in the interior region of our silica film 
(dashed) and in the surface region (solid bold) (see text for 
the definition of the regions). 
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FIG. 3. O-Si-O and Si-O-Si angle distributions in a bulk 
silica sample (circles), the interior (squares) and surface (di- 
amonds) regions of our silica film. 
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FIG. 4. (a) Plot of P„ versus n in a bulk silica sample, in 
the initial surface region (not relaxed) and in a 42 ps annealed 
surface region (relax;ed) where Pn is the probability that an 
atom in a given region is a member of an n-fold ring, 
(b) Plot of < cos^ d > as a function of the ring size n, where 
d is the angle between the direction perpendicular to the sur- 
face and the direction perpendicular to the ring. The dashed 
line represents < cos^ ^ >= i corresponding to an isotropic 
orientation of the rings. 
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FIG. 5. Plot of the mean-squared displacement (MSD) of 

the bridging oxygens (BO) and non-bridging oxygens (NBO) 
in the interior and surface regions as a function of the an- 
nealing time. We have distinguished the parallel component 
of the MSD 1.5 x (x'^ + y^) (dashed and labeled xy) and the 
perpendicular component 3 x {z^) (solid and labeled z). 
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